# Load packages
library(MASS)
library(ggplot2)
library(stringr)
library(tidyr)
library(plyr)
library(gridExtra)
library(dplyr)
library(stargazer)
library(margins)
library(reshape2)
library(zoo)
library(interactions)
library(lme4)
library(readr)
library(xtable)
library(RItools)
library(lmtest)
library(dotwhisker)
library(sciplot)

# The root for the directory for saving PDFs is currently "/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/[filename]".


# clear
rm(list=ls())

setwd("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/")

geo.d <- read.csv(file="~/Dropbox/Complementarity ZC/Georgia Data Sept 2019/geo_working_2021_07_02.csv",head=TRUE,sep=",")   # The setup file for this is: fall19data_setup_2021_09_14.do
phl.d <- read.csv(file="~/Dropbox/Complementarity ZC/Survey_PHL/PHL data/PHL_processed_2021_07_12.csv",head=TRUE,sep=",")	# The setup file for this is: setup_analysis_phl_09_14_21.do
zaf.d <- read.csv(file="~/Dropbox/Complementarity ZC/Survey_ZA/ZAF data/ZAF_processed_2021_07_12.csv",head=TRUE,sep=",")	# The setup file for this is: setup_analysis_zaf_03_02_22.do
usa.d <- read.csv(file="~/Dropbox/Complementarity ZC/Survey_US/Lucid data Mar 2021/us_working_2021_04_26.csv",head=TRUE,sep=",")  	# The setup file for this is: usa_datasetup_analysis_2021_09_14.do`
	usa.d$treatment <- usa.d$tmt_compl
isr.d <- read.csv(file="~/Dropbox/Complementarity ZC/Survey_ISR/ISR data/ISR_processed_2021_08_16.csv",head=TRUE,sep=",")		# The setup file for this is: isr_setup_analysis_09_14_21.do


###
# Creating lists of control variables
###
#		The country specific control lists have two versions, one without anything that could be considered post-treatment (eg legal system) and one without

phl.ctrls <- c("female","age","education_postsec","income_26thto75th","income_above75th","cebuano","tagalog","catholic","pdpl","newshours_morethan4","metromanilancr")
phl.ctrls.other <- c("language_tagalog")
phl.ctrls.other2 <- c("legsyscapable_bin","language_tagalog")
phl.ctrls.psplit <- c("female","age","education_postsec","income_26thto75th","income_above75th","cebuano","tagalog","catholic","newshours_morethan4","metromanilancr")

zaf.ctrls <- c("female","age","education_postsec","income_26thto75th","income_above75th","English","Zulu","Xhosa","christian","party_anc","party_da","newshours_morethan5","gauteng")
zaf.ctrls.other <- c("salegsyscapable_bin")
zaf.ctrls.other2 <- c("salegsyscapable_bin","otherlegsyscapable_bin")

usa.ctrls <- c("female","age","baorhigher","hhi","region_midwest","region_south","region_west","white_bin","demrep_ordered","newshours_morethan6")
usa.ctrls.other <- c("io_sum")
usa.ctrls.other2 <- c("uslegalsystem_bin","dominst_sum","io_sum")

geo.ctrls <- c("female","age","edu_postsecond","enoughmoney_scale","capital","urban")
geo.ctrls.other <- c("displaced")
geo.ctrls.other2 <- c("legalsystem_capable","displaced")


isr.ctrls <- c("female","age","education_postsec","income_26thto75th","income_above75th","leftright","jewish","newshours_morethan4")
isr.ctrls.other <- c("isrlegsyscapable_bin")

###########################
#### Setup for Table 1 ####
###########################

### OLS regressions
### No controls, just regression of support on treatment
m1.geo.ols <-lm(invsupport_agree ~ treatment, data = geo.d)
m1.phl.ols <-lm(iccsupport_bin ~ treatment, data = phl.d)
m1.usa.ols <-lm(suppinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d)
m1.usa.red.ols <-lm(suppinv_bin ~ treatment, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0))
m1.zaf.ols <-lm(iccsupport_bin ~ treatment, data = zaf.d)
m1.isr.ols <-lm(iccsupport_bin ~ treatment, data = isr.d)

summary(m1.geo.ols)
summary(m1.phl.ols)
summary(m1.usa.red.ols)
summary(m1.zaf.ols)
summary(m1.isr.ols)

### With "standard" controls, ones that aren't country-specific
m1.geo.ctrl.ols <-lm(as.formula(paste("invsupport_agree", paste(c("treatment",geo.ctrls), collapse=" + "), sep=" ~ ")), data = geo.d)
m1.phl.ctrl.ols <-lm(as.formula(paste("iccsupport_bin", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d)
m1.usa.ctrl.ols <-lm(as.formula(paste("suppinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = usa.d)
m1.usa.red.ctrl.ols <-lm(as.formula(paste("suppinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0))
m1.zaf.ctrl.ols <-lm(as.formula(paste("iccsupport_bin", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d)
m1.isr.ctrl.ols <-lm(as.formula(paste("iccsupport_bin", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d)

summary(m1.geo.ctrl.ols)
summary(m1.phl.ctrl.ols)
summary(m1.usa.red.ctrl.ols)
summary(m1.zaf.ctrl.ols)
summary(m1.isr.ctrl.ols)


### Logit regressions
### No controls, just regression of support on treatment
m1.geo.log <-glm(invsupport_agree ~ treatment, data = geo.d, family = "binomial")
m1.phl.log <-glm(iccsupport_bin ~ treatment, data = phl.d, family = "binomial")
m1.usa.log <-glm(suppinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d, family = "binomial")
m1.usa.red.log <-glm(suppinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), family = "binomial")
m1.zaf.log <-glm(iccsupport_bin ~ treatment, data = zaf.d, family = "binomial")
m1.isr.log <-glm(iccsupport_bin ~ treatment, data = isr.d, family = "binomial")

summary(m1.geo.log)
summary(m1.phl.log)
summary(m1.usa.red.log)
summary(m1.zaf.log)
summary(m1.isr.log)


### With "standard" controls, ones that aren't country-specific
m1.geo.ctrl.log <-glm(as.formula(paste("invsupport_agree", paste(c("treatment",geo.ctrls), collapse=" + "), sep=" ~ ")), data = geo.d, family = "binomial")
m1.phl.ctrl.log <-glm(as.formula(paste("iccsupport_bin", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d, family = "binomial")
m1.usa.ctrl.log <-glm(as.formula(paste("suppinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = usa.d, family = "binomial")
m1.usa.red.ctrl.log <-glm(as.formula(paste("suppinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), family = "binomial")
m1.zaf.ctrl.log <-glm(as.formula(paste("iccsupport_bin", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d, family = "binomial")
m1.isr.ctrl.log <-glm(as.formula(paste("iccsupport_bin", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d, family = "binomial")

summary(m1.geo.ctrl.log)
summary(m1.phl.ctrl.log)
summary(m1.usa.red.ctrl.log)
summary(m1.zaf.ctrl.log)
summary(m1.isr.ctrl.log)


### Logit regressions
### No controls, just regression of support on treatment
m2.geo.log <-glm(georgiabetter_agree ~ treatment, data = geo.d, family = "binomial")
m2.phl.log <-glm(phlsupport_bin ~ treatment, data = phl.d, family = "binomial")
m2.usa.log <-glm(usinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d, family = "binomial")
m2.usa.red.log <-glm(usinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), family = "binomial")
m2.zaf.log <-glm(countryinvsupport_bin ~ treatment, data = zaf.d, family = "binomial")
m2.isr.log <-glm(countryinvsupport_bin ~ treatment, data = isr.d, family = "binomial")

### With "standard" controls, ones that aren't country-specific
m2.geo.ctrl.log <-glm(as.formula(paste("georgiabetter_agree", paste(c("treatment",geo.ctrls), collapse=" + "), sep=" ~ ")), data = geo.d, family = "binomial")
m2.phl.ctrl.log <-glm(as.formula(paste("phlsupport_bin", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d, family = "binomial")
m2.usa.ctrl.log <-glm(as.formula(paste("usinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = usa.d, family = "binomial")
m2.usa.red.ctrl.log <-glm(as.formula(paste("usinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), family = "binomial")
m2.zaf.ctrl.log <-glm(as.formula(paste("countryinvsupport_bin", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d, family = "binomial")
m2.isr.ctrl.log <-glm(as.formula(paste("countryinvsupport_bin", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d, family = "binomial")

### Ordered Logit regressions
geo.d$iccsupport_num <- factor(geo.d$iccsupport_num)
phl.d$iccsupport_num <- factor(phl.d$iccsupport_num)
usa.d$suppiccinv_num5 <- factor(usa.d$suppiccinv_num5)
zaf.d$iccsupport_num <- factor(zaf.d$iccsupport_num)
isr.d$iccsupport_num <- factor(isr.d$iccsupport_num)

### No controls, just regression of support on treatment
m1.geo.olog <-polr(iccsupport_num ~ treatment, data = geo.d, Hess = TRUE)
m1.phl.olog <-polr(iccsupport_num ~ treatment, data = phl.d, Hess = TRUE)
m1.usa.olog <-polr(suppiccinv_num5 ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d, Hess = TRUE)
m1.usa.red.olog <-polr(suppiccinv_num5 ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), Hess = TRUE)
m1.zaf.olog <-polr(iccsupport_num ~ treatment, data = zaf.d, Hess = TRUE)
m1.isr.olog <-polr(iccsupport_num ~ treatment, data = isr.d, Hess = TRUE)

summary(m1.geo.olog)
summary(m1.phl.olog)
summary(m1.usa.red.olog)
summary(m1.zaf.olog)
summary(m1.isr.olog)


### With "standard" controls, ones that aren't country-specific
m1.geo.ctrl.olog <-polr(as.formula(paste("iccsupport_num", paste(c("treatment",geo.ctrls), collapse=" + "), sep=" ~ ")), data = geo.d, Hess = TRUE)
m1.phl.ctrl.olog <-polr(as.formula(paste("iccsupport_num", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d, Hess = TRUE)
m1.usa.ctrl.olog <-polr(as.formula(paste("suppiccinv_num5", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = usa.d, Hess = TRUE)
m1.usa.red.ctrl.olog <-polr(as.formula(paste("suppiccinv_num5", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), Hess = TRUE)
m1.zaf.ctrl.olog <-polr(as.formula(paste("iccsupport_num", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d, Hess = TRUE)
m1.isr.ctrl.olog <-polr(as.formula(paste("iccsupport_num", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d, Hess = TRUE)

### OLS regressions
### No controls, just regression of support on treatment
###   These are used for the p values in the all-country beta figure
m2.geo.ols <-lm(georgiabetter_agree ~ treatment, data = geo.d)
m2.phl.ols <-lm(phlsupport_bin ~ treatment, data = phl.d)
m2.usa.ols <-lm(usinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d)
m2.usa.red.ols <-lm(usinv_bin ~ treatment, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0))
m2.zaf.ols <-lm(countryinvsupport_bin ~ treatment, data = zaf.d)
m2.isr.ols <-lm(countryinvsupport_bin ~ treatment, data = isr.d)
### With "standard" controls, ones that aren't country-specific
m2.geo.ctrl.ols <-lm(as.formula(paste("georgiabetter_agree", paste(c("treatment",geo.ctrls), collapse=" + "), sep=" ~ ")), data = geo.d)
m2.phl.ctrl.ols <-lm(as.formula(paste("phlsupport_bin", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d)
m2.usa.ctrl.ols <-lm(as.formula(paste("usinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = usa.d)
m2.usa.red.ctrl.ols <-lm(as.formula(paste("usinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0))
m2.zaf.ctrl.ols <-lm(as.formula(paste("countryinvsupport_bin", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d)
m2.isr.ctrl.ols <-lm(as.formula(paste("countryinvsupport_bin", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d)

geo.d$domsupport_num <- factor(geo.d$gabetter_num)
phl.d$domsupport_num <- factor(phl.d$phlsupport_num)
usa.d$domsupport_num5 <- factor(usa.d$suppusinv_num5)
zaf.d$domsupport_num <- factor(zaf.d$countryinvsupport_num)
isr.d$domsupport_num <- factor(isr.d$countryinvsupport_num)

### No controls, just regression of support on treatment
m2.geo.olog <-polr(domsupport_num ~ treatment, data = geo.d, Hess = TRUE)
m2.phl.olog <-polr(domsupport_num ~ treatment, data = phl.d, Hess = TRUE)
m2.usa.olog <-polr(domsupport_num5 ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d, Hess = TRUE)
m2.usa.red.olog <-polr(domsupport_num5 ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), Hess = TRUE)
m2.zaf.olog <-polr(domsupport_num ~ treatment, data = zaf.d, Hess = TRUE)
m2.isr.olog <-polr(domsupport_num ~ treatment, data = isr.d, Hess = TRUE)

### With "standard" controls, ones that aren't country-specific
m2.geo.ctrl.olog <-polr(as.formula(paste("domsupport_num", paste(c("treatment",geo.ctrls), collapse=" + "), sep=" ~ ")), data = geo.d, Hess = TRUE)
m2.phl.ctrl.olog <-polr(as.formula(paste("domsupport_num", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d, Hess = TRUE)
m2.usa.ctrl.olog <-polr(as.formula(paste("domsupport_num5", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = usa.d, Hess = TRUE)
m2.usa.red.ctrl.olog <-polr(as.formula(paste("domsupport_num5", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), Hess = TRUE)
m2.zaf.ctrl.olog <-polr(as.formula(paste("domsupport_num", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d, Hess = TRUE)
m2.isr.ctrl.olog <-polr(as.formula(paste("domsupport_num", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d, Hess = TRUE)


#################################
#################################
#################################
############# START HERE ########
#################################
#################################
#################################


#####
# Some stargazer commands return an error message saying "Error in if (is.na(s)) { : the condition has length > 1".  If this happens, then you need to also run this code.
# h/t Alexey Knorre for this fix.  https://gist.github.com/alexeyknorre/b0780836f4cec04d41a863a683f91b53

## Quick fix for stargazer <= 5.2.3 is.na() issue with long model names in R >= 4.2
  # Unload stargazer if loaded
#detach("package:stargazer",unload=T)
  # Delete it
#remove.packages("stargazer")
  # Download the source
#download.file("https://cran.r-project.org/src/contrib/stargazer_5.2.3.tar.gz", destfile = "stargazer_5.2.3.tar.gz")
  # Unpack
#untar("stargazer_5.2.3.tar.gz")
  # Read the sourcefile with .inside.bracket fun
#stargazer_src <- readLines("stargazer/R/stargazer-internal.R")
  # Move the length check 5 lines up so it precedes is.na(.)
#stargazer_src[1990] <- stargazer_src[1995]
#stargazer_src[1995] <- ""
  # Save back
#writeLines(stargazer_src, con="stargazer/R/stargazer-internal.R")
  # Compile and install the patched package
#install.packages("stargazer", repos = NULL, type="source")
#library(stargazer)



#############
## Table 1 ##
#############

# Regressions for N Values

### GEO ###
# Full regression table, GEO.  Also used for labelling stars in the legend

stargazer(m2.geo.ols, m2.geo.ctrl.ols, m2.geo.log, m2.geo.ctrl.log, m2.geo.olog, m2.geo.ctrl.olog, type = "latex",
          title = "Effect of Treatment on Support for Domestic Investigation (GEO)",
          label="tab:tmtdomsuppgeo",
          omit.stat=c("f", "ser","adj.rsq","aic"))

### USA ###
# Full regression table, USA.  Also used for labelling stars in the legend
stargazer(m1.usa.ols, m1.usa.ctrl.ols, m1.usa.log, m1.usa.ctrl.log, m1.usa.olog, m1.usa.ctrl.olog, type = "latex",
          title = "Effect of Treatment on Support for ICC Investigation (USA; All Treatments)",
          label="tab:tmticcsuppusa",
          omit.stat=c("f", "ser","adj.rsq","aic"))


### PHL ###
# Full regression table, PHL.  Also used for labelling stars in the legend
stargazer(m2.phl.ols, m2.phl.ctrl.ols, m2.phl.log, m2.phl.ctrl.log, m2.phl.olog, m2.phl.ctrl.olog, type = "latex",
          title = "Effect of Treatment on Support for Domestic Investigation (PHL)",
          label="tab:tmtdomsuppphl",
          omit.stat=c("f", "ser","adj.rsq"))

### ZAF ###
# Full regression table, ZAF.  Also used for labelling stars in the legend
stargazer(m1.zaf.ols, m1.zaf.ctrl.ols, m1.zaf.log, m1.zaf.ctrl.log, m1.zaf.olog, m1.zaf.ctrl.olog, type = "latex",
          title = "Effect of Treatment on Support for ICC Investigation (ZAF)",
          label="tab:tmticcsuppzaf",
          omit.stat=c("f", "ser","adj.rsq","aic"))

### ISR ###
# Full regression table, ISR.  Also used for labelling stars in the legend
stargazer(m1.isr.ols, m1.isr.ctrl.ols, m1.isr.log, m1.isr.ctrl.log, m1.isr.olog, m1.isr.ctrl.olog, type = "latex",
          title = "Effect of Treatment on Support for ICC Investigation (ISR)",
          label="tab:tmticcsuppisr",
          omit.stat=c("f", "ser","adj.rsq","aic"))

##############
## Figure 2 ##
##############

### BETA FIGURE, ALL COUNTRIES, Support for ICC
nbyti<- as.matrix(aggregate(iccsupport_bin ~ treatment, data = isr.d, FUN=function(x) sum(!is.na(x)) ))
n1<-nbyti[1,2] # Control, ISR
n2<-nbyti[2,2] # Treatment, ISR
nbyti<- as.matrix(aggregate(invsupport_agree ~ treatment, data = geo.d, FUN=function(x) sum(!is.na(x)) ))
n3<-nbyti[1,2] # Control, GEO
n4<-nbyti[2,2] # Treatment, GEO
nbyti<- as.matrix(aggregate(iccsupport_bin ~ treatment, data = phl.d, FUN=function(x) sum(!is.na(x)) ))
n5<-nbyti[1,2] # Control, PHL
n6<-nbyti[2,2] # Treatment, PHL
nbyti<- as.matrix(aggregate(suppinv_bin ~ treatment, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), FUN=function(x) sum(!is.na(x)) ))
n7<-nbyti[1,2] # Control, USA
n8<-nbyti[2,2] # Treatment, USA
nbyti<- as.matrix(aggregate(iccsupport_bin ~ treatment, data = zaf.d, FUN=function(x) sum(!is.na(x)) ))
n9<-nbyti[1,2] # Control, ZAF
n10<-nbyti[2,2] # Treatment, ZAF

d2 <- isr.d[ which(isr.d$iccsupport_bin =="1"), ]
abyti<- as.matrix(aggregate(iccsupport_bin ~ treatment, data=d2, FUN=function(x) sum(!is.na(x)) ))
a1<-abyti[1,2]	# Approve, Control, ISR
a2<-abyti[2,2]	# Approve, Treatment, ISR
d2 <- geo.d[ which(geo.d$invsupport_agree =="1"), ]
abyti<- as.matrix(aggregate(invsupport_agree ~ treatment, data=d2, FUN=function(x) sum(!is.na(x)) ))
a3<-abyti[1,2]	# Approve, Control, GEO
a4<-abyti[2,2]	# Approve, Treatment, GEO
d2 <- phl.d[ which(phl.d$iccsupport_bin =="1"), ]
abyti<- as.matrix(aggregate(iccsupport_bin ~ treatment, data=d2, FUN=function(x) sum(!is.na(x)) ))
a5<-abyti[1,2]	# Approve, Control, PHL
a6<-abyti[2,2]	# Approve, Treatment, PHL
d2 <- usa.d[ which(usa.d$suppinv_bin =="1" & usa.d$tmt_bias == 0 & usa.d$tmt_biascompl == 0 & usa.d$tmt_biashr == 0), ]
abyti<- as.matrix(aggregate(suppinv_bin ~ treatment, data=d2, FUN=function(x) sum(!is.na(x)) ))
a7<-abyti[1,2]	# Approve, Control, USA
a8<-abyti[2,2]	# Approve, Treatment, USA
d2 <- zaf.d[ which(zaf.d$iccsupport_bin =="1"), ]
abyti<- as.matrix(aggregate(iccsupport_bin ~ treatment, data=d2, FUN=function(x) sum(!is.na(x)) ))
a9<-abyti[1,2]	# Approve, Control, ZAF
a10<-abyti[2,2]	# Approve, Treatment, ZAF

n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
pi_3 <- as.matrix(rbeta(n, a3+0.5, n3-a3+0.5))
pi_4 <- as.matrix(rbeta(n, a4+0.5, n4-a4+0.5))
pi_5 <- as.matrix(rbeta(n, a5+0.5, n5-a5+0.5))
pi_6 <- as.matrix(rbeta(n, a6+0.5, n6-a6+0.5))
pi_7 <- as.matrix(rbeta(n, a7+0.5, n7-a7+0.5))
pi_8 <- as.matrix(rbeta(n, a8+0.5, n8-a8+0.5))
pi_9 <- as.matrix(rbeta(n, a9+0.5, n9-a9+0.5))
pi_10 <- as.matrix(rbeta(n, a10+0.5, n10-a10+0.5))

t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
t3 <- as.matrix(rep(3,n, nrow=n, ncol=1))
t4 <- as.matrix(rep(4,n, nrow=n, ncol=1))
t5 <- as.matrix(rep(5,n, nrow=n, ncol=1))
t6 <- as.matrix(rep(6,n, nrow=n, ncol=1))
t7 <- as.matrix(rep(7,n, nrow=n, ncol=1))
t8 <- as.matrix(rep(8,n, nrow=n, ncol=1))
t9 <- as.matrix(rep(9,n, nrow=n, ncol=1))
t10 <- as.matrix(rep(10,n, nrow=n, ncol=1))

tis <- as.matrix(rbind(t1,t2,t3,t4,t5,t6,t7,t8,t9,t10))
pis <- as.matrix(rbind(pi_1,pi_2,pi_3,pi_4,pi_5,pi_6,pi_7,pi_8,pi_9,pi_10))
betas <- as.data.frame(cbind(tis,pis))
betas$treatmentname = factor(betas$V1, labels = c("ISR Ctrl.","ISR Tmt.","GEO Ctrl.","GEO Tmt.","PHL Ctrl.","PHL Tmt.","USA Ctrl.","USA Tmt.","ZAF Ctrl.","ZAF Tmt."))
lineplot.CI(x.factor = betas$treatmentname, response = betas$V2, type="p",
            data = betas, main = "", xlab = "Treatment Condition, Country", ylab = "% Approving Investigation",
            ylim=c(0.20,1.0),
            ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)),
            cex.axis = 1.2,
            cex.lab = 1.5,
            )

dev.print(pdf, "/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/all_approveinv_beta.pdf")


############################
#### Setup for Figure 3 ####
############################

### OLS regressions
### No controls, just regression of support on treatment
m1.geo.ols <-lm(invsupport_agree ~ treatment, data = geo.d)
m1.phl.ols <-lm(iccsupport_bin ~ treatment, data = phl.d)
m1.usa.ols <-lm(suppinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d)
m1.usa.red.ols <-lm(suppinv_bin ~ treatment, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0))
m1.zaf.ols <-lm(iccsupport_bin ~ treatment, data = zaf.d)
m1.isr.ols <-lm(iccsupport_bin ~ treatment, data = isr.d)

summary(m1.geo.ols)
summary(m1.phl.ols)
summary(m1.usa.ols)
summary(m1.zaf.ols)
summary(m1.isr.ols)

### With "standard" controls, ones that aren't country-specific
m1.geo.ctrl.ols <-lm(as.formula(paste("invsupport_agree", paste(c("treatment",geo.ctrls), collapse=" + "), sep=" ~ ")), data = geo.d)
m1.phl.ctrl.ols <-lm(as.formula(paste("iccsupport_bin", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d)
m1.usa.ctrl.ols <-lm(as.formula(paste("suppinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = usa.d)
m1.usa.red.ctrl.ols <-lm(as.formula(paste("suppinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0))
m1.zaf.ctrl.ols <-lm(as.formula(paste("iccsupport_bin", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d)
m1.isr.ctrl.ols <-lm(as.formula(paste("iccsupport_bin", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d)

summary(m1.geo.ctrl.ols)
summary(m1.phl.ctrl.ols)
summary(m1.usa.ols)
summary(m1.zaf.ols)
summary(m1.isr.ols)


### Logit regressions
### No controls, just regression of support on treatment
m1.geo.log <-glm(invsupport_agree ~ treatment, data = geo.d, family = "binomial")
m1.phl.log <-glm(iccsupport_bin ~ treatment, data = phl.d, family = "binomial")
m1.usa.log <-glm(suppinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d, family = "binomial")
m1.usa.red.log <-glm(suppinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), family = "binomial")
m1.zaf.log <-glm(iccsupport_bin ~ treatment, data = zaf.d, family = "binomial")
m1.isr.log <-glm(iccsupport_bin ~ treatment, data = isr.d, family = "binomial")

### With "standard" controls, ones that aren't country-specific
m1.geo.ctrl.log <-glm(as.formula(paste("invsupport_agree", paste(c("treatment",geo.ctrls), collapse=" + "), sep=" ~ ")), data = geo.d, family = "binomial")
m1.phl.ctrl.log <-glm(as.formula(paste("iccsupport_bin", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d, family = "binomial")
m1.usa.ctrl.log <-glm(as.formula(paste("suppinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = usa.d, family = "binomial")
m1.usa.red.ctrl.log <-glm(as.formula(paste("suppinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), family = "binomial")
m1.zaf.ctrl.log <-glm(as.formula(paste("iccsupport_bin", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d, family = "binomial")
m1.isr.ctrl.log <-glm(as.formula(paste("iccsupport_bin", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d, family = "binomial")

### Ordered Logit regressions
geo.d$iccsupport_num <- factor(geo.d$iccsupport_num)
phl.d$iccsupport_num <- factor(phl.d$iccsupport_num)
usa.d$suppiccinv_num5 <- factor(usa.d$suppiccinv_num5)
zaf.d$iccsupport_num <- factor(zaf.d$iccsupport_num)
isr.d$iccsupport_num <- factor(isr.d$iccsupport_num)
### No controls, just regression of support on treatment
m1.geo.olog <-polr(iccsupport_num ~ treatment, data = geo.d, Hess = TRUE)
m1.phl.olog <-polr(iccsupport_num ~ treatment, data = phl.d, Hess = TRUE)
m1.usa.olog <-polr(suppiccinv_num5 ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d, Hess = TRUE)
m1.usa.red.olog <-polr(suppiccinv_num5 ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), Hess = TRUE)
m1.zaf.olog <-polr(iccsupport_num ~ treatment, data = zaf.d, Hess = TRUE)
m1.isr.olog <-polr(iccsupport_num ~ treatment, data = isr.d, Hess = TRUE)

### With "standard" controls, ones that aren't country-specific
m1.geo.ctrl.olog <-polr(as.formula(paste("iccsupport_num", paste(c("treatment",geo.ctrls), collapse=" + "), sep=" ~ ")), data = geo.d, Hess = TRUE)
m1.phl.ctrl.olog <-polr(as.formula(paste("iccsupport_num", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d, Hess = TRUE)
m1.usa.ctrl.olog <-polr(as.formula(paste("suppiccinv_num5", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = usa.d, Hess = TRUE)
m1.usa.red.ctrl.olog <-polr(as.formula(paste("suppiccinv_num5", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), Hess = TRUE)
m1.zaf.ctrl.olog <-polr(as.formula(paste("iccsupport_num", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d, Hess = TRUE)
m1.isr.ctrl.olog <-polr(as.formula(paste("iccsupport_num", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d, Hess = TRUE)
############################


##############
## Figure 3 ##
##############


### Positive Results ###
########################
### ISR plot for ICC DV, all models
dwplot(list(m1.isr.ols, m1.isr.ctrl.ols, m1.isr.log, m1.isr.ctrl.log, m1.isr.olog, m1.isr.ctrl.olog),
       model_order = c("Model 1", "Model 2", "Model 3", "Model 4", "Model 5", "Model 6"), 
       vline = geom_vline(xintercept = 0), vars_order = c("treatment"), dodge_size = 0.9, dot_args = list(size = 3.5),
	  whisker_args = list(size = 1.0),) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("ISR") +
  theme(
    plot.title = element_text(face = "bold", size = 30),
    legend.position = c(0.007, 0.1),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 20),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("Ord. Logit w/ctrl ***", "Ord. Logit **", "Logit w/ctrl ***", "Logit ***", "OLS w/ctrl ***","OLS ***")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_ISR_dv1.pdf")



### NULL ###
############

### GEO plot for ICC DV, all models
dwplot(list(m1.geo.ols, m1.geo.ctrl.ols, m1.geo.log, m1.geo.ctrl.log, m1.geo.olog, m1.geo.ctrl.olog),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment"), dodge_size = 0.9, dot_args = list(size = 3.5),
	  whisker_args = list(size = 1.0),) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("GEO") +
  theme(
    plot.title = element_text(face = "bold", size = 30),
    legend.position = c(0.007, .6),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 20),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 12),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("Ord. Logit w/ctrl", "Ord. Logit", "Logit w/ctrl", "Logit", "OLS w/ctrl","OLS")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_GEO_dv1.pdf")

### PHL plot for ICC DV, all models
dwplot(list(m1.phl.ols, m1.phl.ctrl.ols, m1.phl.log, m1.phl.ctrl.log, m1.phl.olog, m1.phl.ctrl.olog),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment"), dodge_size = 0.9, dot_args = list(size = 3.5),
	  whisker_args = list(size = 1.0),) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("PHL") +
  theme(
    plot.title = element_text(face = "bold", size = 30),
    legend.position = c(0.007, .65),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 20),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 12),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("Ord. Logit w/ctrl", "Ord. Logit", "Logit w/ctrl", "Logit", "OLS w/ctrl","OLS")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_PHL_dv1.pdf")


### (Weak) Negative Results ###
###############################

### USA plot for ICC DV, all models (REDUCED)
dwplot(list(m1.usa.ols, m1.usa.ctrl.ols, m1.usa.log, m1.usa.ctrl.log, m1.usa.olog, m1.usa.ctrl.olog),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment"), dodge_size = 0.9, dot_args = list(size = 3.5),
	  whisker_args = list(size = 1.0),) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("USA") +
  theme(
    plot.title = element_text(face = "bold", size = 30),
    legend.position = c(0.002, .65),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 20),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 12),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("Ord. Logit w/ctrl", "Ord. Logit", "Logit w/ctrl", "Logit", "OLS w/ctrl","OLS")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_USA_dv1.pdf")

### ZAF plot for ICC DV, all models
dwplot(list(m1.zaf.ols, m1.zaf.ctrl.ols, m1.zaf.log, m1.zaf.ctrl.log, m1.zaf.olog, m1.zaf.ctrl.olog),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment"), dodge_size = 0.9, dot_args = list(size = 3.5),
	  whisker_args = list(size = 1.0),) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("ZAF") +
  theme(
    plot.title = element_text(face = "bold", size = 30),
    legend.position = c(0.007, .65),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 20),
    legend.text = element_text(size = 20),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 12),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("Ord. Logit w/ctrl", "Ord. Logit", "Logit w/ctrl *", "Logit *", "OLS w/ctrl *","OLS *")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_ZAF_dv1.pdf")



##################
#### Figure 4 ####
##################

### BETA FIGURE, ALL COUNTRIES, Support for Domestic Investigation
nbyti<- as.matrix(aggregate(countryinvsupport_bin ~ treatment, data = isr.d, FUN=function(x) sum(!is.na(x)) ))
n1<-nbyti[1,2] # Control, ISR
n2<-nbyti[2,2] # Treatment, ISR
nbyti<- as.matrix(aggregate(georgiabetter_agree ~ treatment, data = geo.d, FUN=function(x) sum(!is.na(x)) ))
n3<-nbyti[1,2] # Control, GEO
n4<-nbyti[2,2] # Treatment, GEO
nbyti<- as.matrix(aggregate(phlsupport_bin ~ treatment, data = phl.d, FUN=function(x) sum(!is.na(x)) ))
n5<-nbyti[1,2] # Control, PHL
n6<-nbyti[2,2] # Treatment, PHL
nbyti<- as.matrix(aggregate(usinv_bin ~ treatment, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), FUN=function(x) sum(!is.na(x)) ))
n7<-nbyti[1,2] # Control, USA
n8<-nbyti[2,2] # Treatment, USA
nbyti<- as.matrix(aggregate(countryinvsupport_bin ~ treatment, data = zaf.d, FUN=function(x) sum(!is.na(x)) ))
n9<-nbyti[1,2] # Control, ZAF
n10<-nbyti[2,2] # Treatment, ZAF

d2 <- isr.d[ which(isr.d$countryinvsupport_bin =="1"), ]
abyti<- as.matrix(aggregate(countryinvsupport_bin ~ treatment, data=d2, FUN=function(x) sum(!is.na(x)) ))
a1<-abyti[1,2]	# Approve, Control, ISR
a2<-abyti[2,2]	# Approve, Treatment, ISR
d2 <- geo.d[ which(geo.d$georgiabetter_agree =="1"), ]
abyti<- as.matrix(aggregate(georgiabetter_agree ~ treatment, data=d2, FUN=function(x) sum(!is.na(x)) ))
a3<-abyti[1,2]	# Approve, Control, GEO
a4<-abyti[2,2]	# Approve, Treatment, GEO
d2 <- phl.d[ which(phl.d$phlsupport_bin =="1"), ]
abyti<- as.matrix(aggregate(phlsupport_bin ~ treatment, data=d2, FUN=function(x) sum(!is.na(x)) ))
a5<-abyti[1,2]	# Approve, Control, PHL
a6<-abyti[2,2]	# Approve, Treatment, PHL
d2 <- usa.d[ which(usa.d$usinv_bin =="1" & usa.d$tmt_bias == 0 & usa.d$tmt_biascompl == 0 & usa.d$tmt_biashr == 0), ]
abyti<- as.matrix(aggregate(usinv_bin ~ treatment, data=d2, FUN=function(x) sum(!is.na(x)) ))
a7<-abyti[1,2]	# Approve, Control, USA
a8<-abyti[2,2]	# Approve, Treatment, USA
d2 <- zaf.d[ which(zaf.d$countryinvsupport_bin =="1"), ]
abyti<- as.matrix(aggregate(countryinvsupport_bin ~ treatment, data=d2, FUN=function(x) sum(!is.na(x)) ))
a9<-abyti[1,2]	# Approve, Control, ZAF
a10<-abyti[2,2]	# Approve, Treatment, ZAF

n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
pi_3 <- as.matrix(rbeta(n, a3+0.5, n3-a3+0.5))
pi_4 <- as.matrix(rbeta(n, a4+0.5, n4-a4+0.5))
pi_5 <- as.matrix(rbeta(n, a5+0.5, n5-a5+0.5))
pi_6 <- as.matrix(rbeta(n, a6+0.5, n6-a6+0.5))
pi_7 <- as.matrix(rbeta(n, a7+0.5, n7-a7+0.5))
pi_8 <- as.matrix(rbeta(n, a8+0.5, n8-a8+0.5))
pi_9 <- as.matrix(rbeta(n, a9+0.5, n9-a9+0.5))
pi_10 <- as.matrix(rbeta(n, a10+0.5, n10-a10+0.5))

t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
t3 <- as.matrix(rep(3,n, nrow=n, ncol=1))
t4 <- as.matrix(rep(4,n, nrow=n, ncol=1))
t5 <- as.matrix(rep(5,n, nrow=n, ncol=1))
t6 <- as.matrix(rep(6,n, nrow=n, ncol=1))
t7 <- as.matrix(rep(7,n, nrow=n, ncol=1))
t8 <- as.matrix(rep(8,n, nrow=n, ncol=1))
t9 <- as.matrix(rep(9,n, nrow=n, ncol=1))
t10 <- as.matrix(rep(10,n, nrow=n, ncol=1))

tis <- as.matrix(rbind(t1,t2,t3,t4,t5,t6,t7,t8,t9,t10))
pis <- as.matrix(rbind(pi_1,pi_2,pi_3,pi_4,pi_5,pi_6,pi_7,pi_8,pi_9,pi_10))
betas <- as.data.frame(cbind(tis,pis))
betas$treatmentname = factor(betas$V1, labels = c("ISR Ctrl.","ISR Tmt.","GEO Ctrl.","GEO Tmt.","PHL Ctrl.","PHL Tmt.","USA Ctrl.","USA Tmt.","ZAF Ctrl.","ZAF Tmt."))
lineplot.CI(x.factor = betas$treatmentname, response = betas$V2, type="p",
            data = betas, main = "", xlab = "Treatment Condition, Country", ylab = "% Approving Dom. Investigation",
            ylim=c(0.30,0.9),
            ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))

dev.print(pdf, "/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/all_approvedomesticinv_beta.pdf")



##################
#### Figure 5 ####
##################

###	Beta Figure: Investigation DV, Full Sample (Excluding DKRTA)
###		m1.usa.ols <-lm(suppinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d)

nbyti<- as.matrix(aggregate(suppinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data=usa.d, FUN=function(x) sum(!is.na(x)) ))
n1<-nbyti[1,5] # Control
n2<-nbyti[2,5] # Treatment
n3<-nbyti[3,5] # Bias
n4<-nbyti[4,5] # Bias + C
n5<-nbyti[5,5] # Bias + HR
d2 <- usa.d[ which(usa.d$suppinv_bin =="1"), ]
abyti<- as.matrix(aggregate(suppinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data=d2, FUN=function(x) sum(!is.na(x)) ))
a1<-abyti[1,5]
a2<-abyti[2,5]
a3<-abyti[3,5]
a4<-abyti[4,5]
a5<-abyti[5,5]
n<-5000
pi_1 <- as.matrix(rbeta(n, a1+0.5, n1-a1+0.5))
pi_2 <- as.matrix(rbeta(n, a2+0.5, n2-a2+0.5))
pi_3 <- as.matrix(rbeta(n, a3+0.5, n3-a3+0.5))
pi_4 <- as.matrix(rbeta(n, a4+0.5, n4-a4+0.5))
pi_5 <- as.matrix(rbeta(n, a5+0.5, n5-a5+0.5))
t1 <- as.matrix(rep(1,n, nrow=n, ncol=1))
t2 <- as.matrix(rep(2,n, nrow=n, ncol=1))
t3 <- as.matrix(rep(3,n, nrow=n, ncol=1))
t4 <- as.matrix(rep(4,n, nrow=n, ncol=1))
t5 <- as.matrix(rep(5,n, nrow=n, ncol=1))
tis <- as.matrix(rbind(t1,t2,t3,t4,t5))
pis <- as.matrix(rbind(pi_1,pi_2,pi_3,pi_4,pi_5))
betas <- as.data.frame(cbind(tis,pis))
betas$treatmentname = factor(betas$V1, labels = c("Control","Compl.","Bias","Bias + Compl.","Bias + HR"))
lineplot.CI(x.factor = betas$treatmentname, response = betas$V2, type="p",
            data = betas, main = "", xlab = "Treatment Group", ylab = "% Approving Investigation",
            ylim=c(0.60,0.8),
            ci.fun = function(x) quantile(x, probs=c(0.05, 0.95)))
dev.print(pdf, "/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/usa_approveinv_beta.pdf")





##############################
##################
#### Appendix B ####
##################
##############################

###################
#### Table A1 ####
###################

### This is just text from the various treatments.

##############################
##################
#### Appendix D ####
##################
##############################



###################
#### Figure A1 ####
###################

#### BALANCE ####
#################

xBalance(as.formula(paste("treatment", paste(c(usa.ctrls), collapse=" + "), sep=" ~ ")), data= subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0),
         report="all", na.rm = TRUE)

xBalance(as.formula(paste("tmt_bias", paste(c(usa.ctrls), collapse=" + "), sep=" ~ ")), data= subset(usa.d, treatment == 0 & tmt_biascompl == 0 & tmt_biashr == 0),
         report="all", na.rm = TRUE)

xBalance(as.formula(paste("tmt_biascompl", paste(c(usa.ctrls), collapse=" + "), sep=" ~ ")), data= subset(usa.d, treatment == 0 & tmt_bias == 0 & tmt_biashr == 0),
         report="all", na.rm = TRUE)

xBalance(as.formula(paste("tmt_biashr", paste(c(usa.ctrls), collapse=" + "), sep=" ~ ")), data= subset(usa.d, treatment == 0 & tmt_bias == 0 & tmt_biascompl == 0),
         report="all", na.rm = TRUE)

### Imbalance is due to the newshour question; underrepresented in the control group

m1.usa.ctrl.ols.news <-lm(as.formula(paste("suppinv_bin", paste(c("treatment*newshours_morethan6","tmt_bias*newshours_morethan6","tmt_biascompl*newshours_morethan6","tmt_biashr*newshours_morethan6",usa.ctrls), collapse=" + "), sep=" ~ ")), data = usa.d)

stargazer(m1.usa.ctrl.ols.news , type = "latex",
          label="tab:newsint",
          omit.stat=c("f", "ser","adj.rsq"))

usa.balance <-RItools::xBalance(
  fmla=as.formula(paste("treatment", paste(c(usa.ctrls), collapse=" + "), sep=" ~ ")),
  na.rm=TRUE,
  data= subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0),
  report = c("std.diffs","p.values","chisquare.test"))

usa.balance.plot <-usa.balance $results %>% 
  data.frame() %>% 
  tibble::rownames_to_column() %>% 
  rename(std.diff=std.diff.unstrat,p=p.unstrat) %>% 
  mutate(rowname=stringr::str_replace(rowname,"female","Female ")) %>%
  mutate(rowname=stringr::str_replace(rowname,"age","Age ")) %>%
  mutate(rowname=stringr::str_replace(rowname,"baorhigher","BA or higher")) %>%
  mutate(rowname=stringr::str_replace(rowname,"hhi","Household Income")) %>%
  mutate(rowname=stringr::str_replace(rowname,"region_midwest","Midwest")) %>%
  mutate(rowname=stringr::str_replace(rowname,"region_south","South")) %>%
  mutate(rowname=stringr::str_replace(rowname,"region_west","West")) %>%
  mutate(rowname=stringr::str_replace(rowname,"white_bin","White")) %>%
  mutate(rowname=stringr::str_replace(rowname,"demrep_ordered","Dem Rep Scale ")) %>%
  mutate(rowname=stringr::str_replace(rowname,"newshours_morethan6","News more than 6 hours")) %>%
  ggplot(aes(std.diff,rowname,color=interaction(p<0.05,p<0.1)))+
  geom_vline(xintercept=0,linetype="dashed",color="gray25")+
  geom_point()+
  scale_color_viridis_d(labels=c("p>0.10","p<0.1","p<0.05"),begin = 0.2,end=0.8)+
  coord_cartesian(xlim = c(-0.20,0.20))+
  theme_bw()+
  theme(legend.position="bottom",axis.text.y=element_text(size=5))+
  labs(x="Standardized Difference",
       y="Variable",
       color="Significance",
       title="Treatment (Control or Complementarity)")
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/balancetable_usa.pdf")
usa.balance.plot

#### Sensitivity Analysis #####
###############################

# USA Sensitivity testing
library(sensemakr)

# Regress outcome on treatment and observables
m1.usa.red.ctrl.ols <-lm(as.formula(paste("suppinv_bin", paste(c("treatment",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0))

# Regress treatment on observables
model.dx <- lm(treatment ~ female+age+baorhigher+hhi+region_midwest+region_south+region_west+white_bin+demrep_ordered+newshours_morethan6, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0))

# Calculate partial R2s
r2yx.d <- partial_r2(m1.usa.red.ctrl.ols, covariates = "newshours_morethan6")
r2dx   <- partial_r2(model.dx, covariates = "newshours_morethan6")

# Informal adjustment
informal_adjusted_estimate.anc <- adjusted_estimate(m1.usa.red.ctrl.ols, 
                                                    treatment = "treatment", 
                                                    r2dz.x = r2dx, 
                                                    r2yz.dx = r2yx.d)

# Formal bounds
formal_bound <- ovb_bounds(model = m1.usa.red.ctrl.ols, 
                           treatment = "treatment", 
                           benchmark_covariates = "newshours_morethan6", 
                           kd = 1, ky = 1)                                                

ovb_contour_plot(m1.usa.red.ctrl.ols,  
                 treatment = "treatment",
                 lim = .01)

# adds proper bound
add_bound_to_contour(bounds = formal_bound,
                     bound_label = "News Hours >= 6")
dev.print(pdf, "/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/contours_usa.pdf")




###################
#### Figure A2 ####
###################

#### BALANCE ####
#################

xBalance(as.formula(paste("treatment", paste(c(zaf.ctrls), collapse=" + "), sep=" ~ ")), data=zaf.d,
         report="all", na.rm = TRUE)

zaf.balance <-RItools::xBalance(
  fmla=as.formula(paste("treatment", paste(c(zaf.ctrls), collapse=" + "), sep=" ~ ")),
  na.rm=TRUE,
  data= zaf.d,
  report = c("std.diffs","p.values","chisquare.test"))

zaf.balance.plot <-zaf.balance $results %>% 
  data.frame() %>% 
  tibble::rownames_to_column() %>% 
  rename(std.diff=std.diff.unstrat,p=p.unstrat) %>% 
  mutate(rowname=stringr::str_replace(rowname,"female","Female ")) %>%
  mutate(rowname=stringr::str_replace(rowname,"age","Age ")) %>%
  mutate(rowname=stringr::str_replace(rowname,"education_postsec","Post-Sec. Ed.")) %>%
  mutate(rowname=stringr::str_replace(rowname,"income_above75th","Income above 75th")) %>%
  mutate(rowname=stringr::str_replace(rowname,"income_26thto75th","Income 26th to 75th")) %>%
  mutate(rowname=stringr::str_replace(rowname,"christian","Christian")) %>%
  mutate(rowname=stringr::str_replace(rowname,"party_anc","ANC")) %>%
  mutate(rowname=stringr::str_replace(rowname,"party_da","Dem. Alliance")) %>%
  mutate(rowname=stringr::str_replace(rowname,"gauteng","Gauteng")) %>%
  mutate(rowname=stringr::str_replace(rowname,"newshours_morethan5","News more than 5 hours")) %>%
  ggplot(aes(std.diff,rowname,color=interaction(p<0.05,p<0.1)))+
  geom_vline(xintercept=0,linetype="dashed",color="gray25")+
  geom_point()+
  scale_color_viridis_d(labels=c("p>0.10","p<0.1","p<0.05"),begin = 0.2,end=0.8)+
  coord_cartesian(xlim = c(-0.10,0.10))+
  theme_bw()+
  theme(legend.position="bottom",axis.text.y=element_text(size=5))+
  labs(x="Standardized Difference",
       y="Variable",
       color="Significance",
       title="Treatment (Control or Complementarity)")
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/balancetable_zaf.pdf")

zaf.balance.plot


#### Sensitivity Analysis #####
###############################

# Formal bounds
# Only doing these for ANC, since it has the highest partial R2 with treatment and just barely the second highest partial R2 for outcome.
formal_bound <- ovb_bounds(model = m1.zaf.ctrl.ols, 
                           treatment = "treatment", 
                           benchmark_covariates = "party_anc", 
                           kd = 1, ky = 1)                                                

ovb_contour_plot(m1.zaf.ctrl.ols,  
                 treatment = "treatment",
                 lim = .005)

# adds proper bound
add_bound_to_contour(bounds = formal_bound,
                     bound_label = "ANC")
dev.print(pdf, "/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/contours_zaf")





##############################
##################
#### Appendix E ####
##################
##############################


##################
#### Table A2 ####
##################

#### Israel ####
summary(m1.isr.ols)
summary(m1.isr.ctrl.ols)
summary(m1.isr.log)
summary(m1.isr.ctrl.log)
summary(m1.isr.olog)
summary(m1.isr.ctrl.olog)

### Georgia ###
summary(m1.geo.ols)
summary(m1.geo.ctrl.ols)
summary(m1.geo.log)
summary(m1.geo.ctrl.log)
summary(m1.geo.olog)
summary(m1.geo.ctrl.olog)

### Philippines ####
summary(m1.phl.ols)
summary(m1.phl.ctrl.ols)
summary(m1.phl.log)
summary(m1.phl.ctrl.log)
summary(m1.phl.olog)
summary(m1.phl.ctrl.olog)

### USA ###
summary(m1.usa.ols)
summary(m1.usa.ctrl.ols)
summary(m1.usa.log)
summary(m1.usa.ctrl.log)
summary(m1.usa.olog)
summary(m1.usa.ctrl.olog)

### South Africa ###
summary(m1.zaf.ols)
summary(m1.zaf.ctrl.ols)
summary(m1.zaf.log)
summary(m1.zaf.ctrl.log)
summary(m1.zaf.olog)
summary(m1.zaf.ctrl.olog)


#############################
#### Setup for Figure A3 ####
#############################
### Logit regressions
### No controls, just regression of support on treatment
m2.geo.log <-glm(georgiabetter_agree ~ treatment, data = geo.d, family = "binomial")
m2.phl.log <-glm(phlsupport_bin ~ treatment, data = phl.d, family = "binomial")
m2.usa.log <-glm(usinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d, family = "binomial")
m2.usa.red.log <-glm(usinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), family = "binomial")
m2.zaf.log <-glm(countryinvsupport_bin ~ treatment, data = zaf.d, family = "binomial")
m2.isr.log <-glm(countryinvsupport_bin ~ treatment, data = isr.d, family = "binomial")

### With "standard" controls, ones that aren't country-specific
m2.geo.ctrl.log <-glm(as.formula(paste("georgiabetter_agree", paste(c("treatment",geo.ctrls), collapse=" + "), sep=" ~ ")), data = geo.d, family = "binomial")
m2.phl.ctrl.log <-glm(as.formula(paste("phlsupport_bin", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d, family = "binomial")
m2.usa.ctrl.log <-glm(as.formula(paste("usinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = usa.d, family = "binomial")
m2.usa.red.ctrl.log <-glm(as.formula(paste("usinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), family = "binomial")
m2.zaf.ctrl.log <-glm(as.formula(paste("countryinvsupport_bin", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d, family = "binomial")
m2.isr.ctrl.log <-glm(as.formula(paste("countryinvsupport_bin", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d, family = "binomial")

### OLS regressions
### No controls, just regression of support on treatment
###   These are used for the p values in the all-country beta figure
m2.geo.ols <-lm(georgiabetter_agree ~ treatment, data = geo.d)
m2.phl.ols <-lm(phlsupport_bin ~ treatment, data = phl.d)
m2.usa.ols <-lm(usinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d)
m2.usa.red.ols <-lm(usinv_bin ~ treatment, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0))
m2.zaf.ols <-lm(countryinvsupport_bin ~ treatment, data = zaf.d)
m2.isr.ols <-lm(countryinvsupport_bin ~ treatment, data = isr.d)

### With "standard" controls, ones that aren't country-specific
m2.geo.ctrl.ols <-lm(as.formula(paste("georgiabetter_agree", paste(c("treatment",geo.ctrls), collapse=" + "), sep=" ~ ")), data = geo.d)
m2.phl.ctrl.ols <-lm(as.formula(paste("phlsupport_bin", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d)
m2.usa.ctrl.ols <-lm(as.formula(paste("usinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = usa.d)
m2.usa.red.ctrl.ols <-lm(as.formula(paste("usinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0))
m2.zaf.ctrl.ols <-lm(as.formula(paste("countryinvsupport_bin", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d)
m2.isr.ctrl.ols <-lm(as.formula(paste("countryinvsupport_bin", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d)

geo.d$domsupport_num <- factor(geo.d$gabetter_num)
phl.d$domsupport_num <- factor(phl.d$phlsupport_num)
usa.d$domsupport_num5 <- factor(usa.d$suppusinv_num5)
zaf.d$domsupport_num <- factor(zaf.d$countryinvsupport_num)
isr.d$domsupport_num <- factor(isr.d$countryinvsupport_num)

### No controls, just regression of support on treatment
m2.geo.olog <-polr(domsupport_num ~ treatment, data = geo.d, Hess = TRUE)
m2.phl.olog <-polr(domsupport_num ~ treatment, data = phl.d, Hess = TRUE)
m2.usa.olog <-polr(domsupport_num5 ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d, Hess = TRUE)
m2.usa.red.olog <-polr(domsupport_num5 ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), Hess = TRUE)
m2.zaf.olog <-polr(domsupport_num ~ treatment, data = zaf.d, Hess = TRUE)
m2.isr.olog <-polr(domsupport_num ~ treatment, data = isr.d, Hess = TRUE)

### With "standard" controls, ones that aren't country-specific
m2.geo.ctrl.olog <-polr(as.formula(paste("domsupport_num", paste(c("treatment",geo.ctrls), collapse=" + "), sep=" ~ ")), data = geo.d, Hess = TRUE)
m2.phl.ctrl.olog <-polr(as.formula(paste("domsupport_num", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d, Hess = TRUE)
m2.usa.ctrl.olog <-polr(as.formula(paste("domsupport_num5", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = usa.d, Hess = TRUE)
m2.usa.red.ctrl.olog <-polr(as.formula(paste("domsupport_num5", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0), Hess = TRUE)
m2.zaf.ctrl.olog <-polr(as.formula(paste("domsupport_num", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d, Hess = TRUE)
m2.isr.ctrl.olog <-polr(as.formula(paste("domsupport_num", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d, Hess = TRUE)
#############################

##################
#### Figure A3 ####
##################


#### Positive Results ####
##########################

### ISR plot for ICC DV, all models
dwplot(list(m2.isr.ols, m2.isr.ctrl.ols, m2.isr.log, m2.isr.ctrl.log, m2.isr.olog, m2.isr.ctrl.olog),
       model_order = c("Model 1", "Model 2", "Model 3", "Model 4", "Model 5", "Model 6"), 
       vline = geom_vline(xintercept = 0), vars_order = c("treatment")) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("ISR") +
  theme(
    plot.title = element_text(face = "bold", size = 20),
    legend.position = c(0.007, .3),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 13),
    legend.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 12),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("Ord. Logit w/ctrl ***", "Ord. Logit **", "Logit w/ctrl ***", "Logit ***", "OLS w/ctrl ***","OLS ***")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_ISR_dv2.pdf")


### ZAF plot for ICC DV, all models
dwplot(list(m2.zaf.ols, m2.zaf.ctrl.ols, m2.zaf.log, m2.zaf.ctrl.log, m2.zaf.olog, m2.zaf.ctrl.olog),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment"),) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("ZAF") +
  theme(
    plot.title = element_text(face = "bold", size = 20),
    legend.position = c(0.007, .3),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 13),
    legend.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 12),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("Ord. Logit w/ctrl ***", "Ord. Logit **", "Logit w/ctrl ***", "Logit ***", "OLS w/ctrl ***","OLS ***")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_ZAF_dv2.pdf")


#### Null ####
##############

### GEO plot for ICC DV, all models
dwplot(list(m2.geo.ols, m2.geo.ctrl.ols, m2.geo.log, m2.geo.ctrl.log, m2.geo.olog, m2.geo.ctrl.olog),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment"),) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("GEO") +
  theme(
    plot.title = element_text(face = "bold", size = 20),
    legend.position = c(0.007, .3),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 13),
    legend.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 12),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("Ord. Logit w/ctrl ***", "Ord. Logit **", "Logit w/ctrl ***", "Logit ***", "OLS w/ctrl ***","OLS ***")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_GEO_dv2.pdf")

### USA plot for ICC DV, all models
dwplot(list(m2.usa.red.ols, m2.usa.red.ctrl.ols, m2.usa.red.log, m2.usa.red.ctrl.log, m2.usa.red.olog, m2.usa.red.ctrl.olog),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment"),) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("USA") +
  theme(
    plot.title = element_text(face = "bold", size = 20),
    legend.position = c(0.007, .3),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 13),
    legend.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 12),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("Ord. Logit w/ctrl ***", "Ord. Logit **", "Logit w/ctrl ***", "Logit ***", "OLS w/ctrl ***","OLS ***")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_USA_dv2.pdf")


#### (Weak) Negative Results ####
#################################

### PHL plot for ICC DV, all models
dwplot(list(m2.phl.ols, m2.phl.ctrl.ols, m2.phl.log, m2.phl.ctrl.log, m2.phl.olog, m2.phl.ctrl.olog),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment"),) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("PHL") +
  theme(
    plot.title = element_text(face = "bold", size = 20),
    legend.position = c(0.007, .3),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 13),
    legend.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 12),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("Ord. Logit w/ctrl ***", "Ord. Logit **", "Logit w/ctrl ***", "Logit ***", "OLS w/ctrl ***","OLS ***")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_PHL_dv2.pdf")

#######

# Tables A3 - A7

stargazer(m2.isr.ols, m2.isr.ctrl.ols, m2.isr.log, m2.isr.ctrl.log, m2.isr.olog, m2.isr.ctrl.olog, type = "latex",
          title = "Effect of Treatment on Support for Domestic Investigation (ISR)",
          label="tab:tmtdomsuppisr",
          omit.stat=c("f", "ser","adj.rsq"))

stargazer(m2.zaf.ols, m2.zaf.ctrl.ols, m2.zaf.log, m2.zaf.ctrl.log, m2.zaf.olog, m2.zaf.ctrl.olog, type = "latex",
          title = "Effect of Treatment on Support for Domestic Investigation (ZAF)",
          label="tab:tmtdomsuppzaf",
          omit.stat=c("f", "ser","adj.rsq"))

stargazer(m2.geo.ols, m2.geo.ctrl.ols, m2.geo.log, m2.geo.ctrl.log, m2.geo.olog, m2.geo.ctrl.olog, type = "latex",
          title = "Effect of Treatment on Support for Domestic Investigation (GEO)",
          label="tab:tmtdomsuppgeo",
          omit.stat=c("f", "ser","adj.rsq"))


stargazer(m2.usa.ols, m2.usa.ctrl.ols, m2.usa.log, m2.usa.ctrl.log, m2.usa.olog, m2.usa.ctrl.olog, type = "latex",
          title = "Effect of Treatment on Support for Domestic Investigation (USA)",
          label="tab:tmtdomsuppusa",
          omit.stat=c("f", "ser","adj.rsq"))


stargazer(m2.phl.ols, m2.phl.ctrl.ols, m2.phl.log, m2.phl.ctrl.log, m2.phl.olog, m2.phl.ctrl.olog, type = "latex",
          title = "Effect of Treatment on Support for Domestic Investigation (PHL)",
          label="tab:tmtdomsuppphl",
          omit.stat=c("f", "ser","adj.rsq"))




##############################
##################
#### Appendix F ####
##################
##############################

###################
#### Setup Figure A4 ####
###################

#####################
# DV 1
#####################
### OLS No controls, just regression of support on treatment
m1.geo.ols.nodk <-lm(invsupport_agree_nodkrta ~ treatment, data = geo.d)
m1.phl.ols.nodk <-lm(iccsupport_bin_nn ~ treatment, data = phl.d)
m1.usa.ols.nodk <-lm(suppinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = subset(usa.d, suppiccinv_num5 != 3))
m1.usa.red.ols.nodk <-lm(suppinv_bin ~ treatment, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0 & suppiccinv_num5 != 3))
m1.zaf.ols.nodk <-lm(iccsupport_bin ~ treatment, data = subset(zaf.d, iccsupport_num != 2))
m1.isr.ols.nodk <-lm(iccsupport_bin_nn ~ treatment, data = isr.d)

### OLS With "standard" controls, ones that aren't country-specific
m1.geo.ctrl.ols.nodk <-lm(as.formula(paste("invsupport_agree_nodkrta", paste(c("treatment",geo.ctrls), collapse=" + "), sep=" ~ ")), data = geo.d)
m1.phl.ctrl.ols.nodk <-lm(as.formula(paste("iccsupport_bin_nn", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d)
m1.usa.ctrl.ols.nodk <-lm(as.formula(paste("suppinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, suppiccinv_num5 != 3))
m1.usa.red.ctrl.ols.nodk <-lm(as.formula(paste("suppinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0 & suppiccinv_num5 != 3))
m1.zaf.ctrl.ols.nodk <-lm(as.formula(paste("iccsupport_bin", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = subset(zaf.d, iccsupport_num != 2))
m1.isr.ctrl.ols.nodk <-lm(as.formula(paste("iccsupport_bin_nn", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d)

### Logit No controls, just regression of support on treatment
m1.geo.log.nodk <-glm(invsupport_agree_nodkrta ~ treatment, data = geo.d, family = "binomial")
m1.phl.log.nodk <-glm(iccsupport_bin_nn ~ treatment, data = phl.d, family = "binomial")
m1.usa.log.nodk <-glm(suppinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = subset(usa.d, suppiccinv_num5 != 3), family = "binomial")
m1.usa.red.log.nodk <-glm(suppinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0 & suppiccinv_num5 != 3), family = "binomial")
m1.zaf.log.nodk <-glm(iccsupport_bin ~ treatment, data = subset(zaf.d, iccsupport_num != 2), family = "binomial")
m1.isr.log.nodk <-glm(iccsupport_bin_nn ~ treatment, data = isr.d, family = "binomial")

### Logit With "standard" controls, ones that aren't country-specific
m1.geo.ctrl.log.nodk <-glm(as.formula(paste("invsupport_agree_nodkrta", paste(c("treatment",geo.ctrls), collapse=" + "), sep=" ~ ")), data = geo.d, family = "binomial")
m1.phl.ctrl.log.nodk <-glm(as.formula(paste("iccsupport_bin_nn", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d, family = "binomial")
m1.usa.ctrl.log.nodk <-glm(as.formula(paste("suppinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, suppiccinv_num5 != 3), family = "binomial")
m1.usa.red.ctrl.log.nodk <-glm(as.formula(paste("suppinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0 & suppiccinv_num5 != 3), family = "binomial")
m1.zaf.ctrl.log.nodk <-glm(as.formula(paste("iccsupport_bin", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = subset(zaf.d, iccsupport_num != 2), family = "binomial")
m1.isr.ctrl.log.nodk <-glm(as.formula(paste("iccsupport_bin_nn", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d, family = "binomial")




###################
#### Figure A4 ####
###################

### Positive ###
################

### ISR plot for ICC DV, all models
dwplot(list(m1.isr.ols.nodk, m1.isr.ctrl.ols.nodk, m1.isr.log.nodk, m1.isr.ctrl.log.nodk),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment")) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("ISR") +
  theme(
    plot.title = element_text(face = "bold", size = 20),
    legend.position = c(0.007, .3),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 13),
    legend.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 12),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("OLS ***","OLS w/ctrl ***","Logit ***","Logit w/ctrl ***")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_ISR_dv1_nodk.pdf")


### NULL ###
############

### GEO plot for ICC DV, all models
dwplot(list(m1.geo.ols.nodk, m1.geo.ctrl.ols.nodk, m1.geo.log.nodk, m1.geo.ctrl.log.nodk),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment"),) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("GEO") +
  theme(
    plot.title = element_text(face = "bold", size = 20),
    legend.position = c(0.007, .3),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 13),
    legend.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 12),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("OLS","OLS w/ctrl","Logit","Logit w/ctrl")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_GEO_dv1_nodk.pdf")

### PHL plot for ICC DV, all models
dwplot(list(m1.phl.ols.nodk, m1.phl.ctrl.ols.nodk, m1.phl.log.nodk, m1.phl.ctrl.log.nodk),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment"),) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("PHL") +
  theme(
    plot.title = element_text(face = "bold", size = 20),
    legend.position = c(0.007, .3),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 13),
    legend.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 12),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("OLS","OLS w/ctrl","Logit","Logit w/ctrl")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_PHL_dv1_nodk.pdf")


### Weak (Negative) ###
#######################

### USA plot for ICC DV, all models (REDUCED)
dwplot(list(m1.usa.red.ols.nodk, m1.usa.red.ctrl.ols.nodk, m1.usa.red.log.nodk, m1.usa.red.ctrl.log.nodk),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment"),) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("USA") +
  theme(
    plot.title = element_text(face = "bold", size = 20),
    legend.position = c(0.007, .3),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 13),
    legend.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 12),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("OLS","OLS w/ctrl","Logit","Logit w/ctrl")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_USA_dv1_nodk.pdf")

### ZAF plot for ICC DV, all models
dwplot(list(m1.zaf.ols.nodk, m1.zaf.ctrl.ols.nodk, m1.zaf.log.nodk, m1.zaf.ctrl.log.nodk),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment"),) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.8,0.1) +
  ggtitle("ZAF") +
  theme(
    plot.title = element_text(face = "bold", size = 20),
    legend.position = c(0.007, .3),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 13),
    legend.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 12),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("OLS *","OLS w/ctrl *","Logit *","Logit w/ctrl *")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_ZAF_dv1_nodk.pdf")


###################
#### Setup for Figure A5 ####
###################

### OLS regressions
### No controls, just regression of support on treatment
m2.geo.ols.nodk <-lm(georgiabetter_agree_nodkrta ~ treatment, data = geo.d)
m2.phl.ols.nodk <-lm(phlsupport_bin_nn ~ treatment, data = phl.d)
m2.usa.ols.nodk <-lm(usinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = subset(usa.d, suppusinv_num5 != 3))
m2.usa.red.ols.nodk <-lm(usinv_bin ~ treatment, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0 & suppusinv_num5 != 3))
m2.zaf.ols.nodk <-lm(countryinvsupport_bin_nn ~ treatment, data = zaf.d)
m2.isr.ols.nodk <-lm(countryinvsupport_bin_nn ~ treatment, data = isr.d)

### With "standard" controls, ones that aren't country-specific
m2.geo.ctrl.ols.nodk <-lm(as.formula(paste("georgiabetter_agree_nodkrta", paste(c("treatment",geo.ctrls), collapse=" + "), sep=" ~ ")), data = geo.d)
m2.phl.ctrl.ols.nodk <-lm(as.formula(paste("phlsupport_bin_nn", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d)
m2.usa.ctrl.ols.nodk <-lm(as.formula(paste("usinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, suppusinv_num5 != 3))
m2.usa.red.ctrl.ols.nodk <-lm(as.formula(paste("usinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0 & suppusinv_num5 != 3))
m2.zaf.ctrl.ols.nodk <-lm(as.formula(paste("countryinvsupport_bin_nn", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d)
m2.isr.ctrl.ols.nodk <-lm(as.formula(paste("countryinvsupport_bin_nn", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d)

### Logit regressions
### No controls, just regression of support on treatment
m2.geo.log.nodk <-glm(georgiabetter_agree_nodkrta ~ treatment, data = geo.d, family = "binomial")
m2.phl.log.nodk <-glm(phlsupport_bin_nn ~ treatment, data = phl.d, family = "binomial")
m2.usa.log.nodk <-glm(usinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = subset(usa.d, suppusinv_num5 != 3), family = "binomial")
m2.usa.red.log.nodk <-glm(usinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0 & suppusinv_num5 != 3), family = "binomial")
m2.zaf.log.nodk <-glm(countryinvsupport_bin_nn ~ treatment, data = zaf.d, family = "binomial")
m2.isr.log.nodk <-glm(countryinvsupport_bin_nn ~ treatment, data = isr.d, family = "binomial")

### With "standard" controls, ones that aren't country-specific
m2.geo.ctrl.log.nodk <-glm(as.formula(paste("georgiabetter_agree_nodkrta", paste(c("treatment",geo.ctrls), collapse=" + "), sep=" ~ ")), data = geo.d, family = "binomial")
m2.phl.ctrl.log.nodk <-glm(as.formula(paste("phlsupport_bin_nn", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d, family = "binomial")
m2.usa.ctrl.log.nodk <-glm(as.formula(paste("usinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, suppusinv_num5 != 3), family = "binomial")
m2.usa.red.ctrl.log.nodk <-glm(as.formula(paste("usinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0 & suppusinv_num5 != 3), family = "binomial")
m2.zaf.ctrl.log.nodk <-glm(as.formula(paste("countryinvsupport_bin_nn", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d, family = "binomial")
m2.isr.ctrl.log.nodk <-glm(as.formula(paste("countryinvsupport_bin_nn", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d, family = "binomial")



###################
#### Figure A5 ####
###################

################
### Positive ###
################
### ISR plot for ICC DV, all models
dwplot(list(m2.isr.ols.nodk, m2.isr.ctrl.ols.nodk, m2.isr.log.nodk, m2.isr.ctrl.log.nodk),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment")) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("ISR") +
  theme(
    plot.title = element_text(face = "bold", size = 20),
    legend.position = c(0.007, 0.3),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 13),
    legend.text = element_text(size = 10),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("OLS ***","OLS w/ctrl ***","Logit ***","Logit w/ctrl ***")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_ISR_dv2_nodk.pdf")


### ZAF plot for ICC DV, all models
dwplot(list(m2.zaf.ols.nodk, m2.zaf.ctrl.ols.nodk, m2.zaf.log.nodk, m2.zaf.ctrl.log.nodk),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment"),) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("ZAF") +
  theme(
    plot.title = element_text(face = "bold", size = 20),
    legend.position = c(0.007, 0.3),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 13),
    legend.text = element_text(size = 10),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("OLS","OLS w/ctrl","Logit","Logit w/ctrl")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_ZAF_dv2_nodk.pdf")


############
### NULL ###
############

### GEO plot for ICC DV, all models
dwplot(list(m2.geo.ols.nodk, m2.geo.ctrl.ols.nodk, m2.geo.log.nodk, m2.geo.ctrl.log.nodk),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment"),) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("GEO") +
  theme(
    plot.title = element_text(face = "bold", size = 20),
    legend.position = c(0.007, 0.3),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 13),
    legend.text = element_text(size = 10),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("OLS","OLS w/ctrl","Logit","Logit w/ctrl")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_GEO_dv2_nodk.pdf")

### USA plot for ICC DV, all models
dwplot(list(m2.usa.red.ols.nodk, m2.usa.red.ctrl.ols.nodk, m2.usa.red.log.nodk, m2.usa.red.ctrl.log.nodk),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment"),) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("USA") +
  theme(
    plot.title = element_text(face = "bold", size = 20),
    legend.position = c(0.007, 0.3),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 13),
    legend.text = element_text(size = 10),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("OLS","OLS w/ctrl","Logit","Logit w/ctrl")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_USA_dv2_nodk.pdf")

################
### Negative ###
################

### PHL plot for ICC DV, all models
dwplot(list(m2.phl.ols.nodk, m2.phl.ctrl.ols.nodk, m2.phl.log.nodk, m2.phl.ctrl.log.nodk),
       vline = geom_vline(xintercept = 0), vars_order = c("treatment"),) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("PHL") +
  theme(
    plot.title = element_text(face = "bold", size = 20),
    legend.position = c(0.007, 0.3),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title = element_text(size = 13),
    legend.text = element_text(size = 10),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 14),
    legend.title.align = .5,
    axis.text.y = element_text(size = 0)
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("OLS ***","OLS w/ctrl**","Logit ***","Logit w/ctrl***")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/Replication Archive/coeffs_PHL_dv2_nodk.pdf")




##############################
##################
#### Appendix G ####
##################
##############################

###################
#### Table A8 #####
###################

man1.phl <-lm(manip1_correct ~ treatment, data = phl.d)
man2.phl <-lm(manip2_correct ~ treatment, data = phl.d)
man1.zaf <-lm(manip1_correct ~ treatment, data = zaf.d)
man2.zaf <-lm(manip2_correct ~ treatment, data = zaf.d)
man1.isr <-lm(manip1_correct ~ treatment, data = isr.d)
man2.isr <-lm(manip2_correct ~ treatment, data = isr.d)

man3.usa <-lm(ccheck1_correct ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d)
man4.usa <-lm(ccheck2_correct ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d)

stargazer(man1.phl, man2.phl, man1.zaf, man2.zaf, man3.usa, man4.usa, man1.isr, man2.isr, type = "latex",
          title = "Effect of Treatment on Manipulation Checks",
          label="tab:mancheck",
          dep.var.labels=c("Q1 PHL","Q2 PHL","Q1 ZAF","Q2 ZAF","Q1 USA","Q2 USA","Q1 ISR","Q2 ISR"), omit.stat=c("f","ser","adj.rsq","aic","ll"))





###################
#### Figure A6 #####
###################

#### Note that these are the correct regressions and dot-whisker plots.  But the order is reversed from the manuscript figure.

usa.d$suppiccinv_num5 <- factor(usa.d$suppiccinv_num5)
usa.d$domsupport_num5 <- factor(usa.d$suppusinv_num5)

usa.d.comp1and2 <- subset(usa.d, complier1and2 == 1)

m1.usa.ols.comp1and2 <-lm(suppinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d.comp1and2)
m1.usa.ctrl.ols.comp1and2 <-lm(as.formula(paste("suppinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = usa.d.comp1and2)
m1.usa.log.comp1and2 <-glm(suppinv_bin ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d.comp1and2, family = "binomial")
m1.usa.ctrl.log.comp1and2 <-glm(as.formula(paste("suppinv_bin", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = usa.d.comp1and2, family = "binomial")
m1.usa.olog.comp1and2 <-polr(suppiccinv_num5 ~ treatment + tmt_bias + tmt_biascompl + tmt_biashr, data = usa.d.comp1and2, Hess = TRUE)
m1.usa.ctrl.olog.comp1and2 <-polr(as.formula(paste("suppiccinv_num5", paste(c("treatment","tmt_bias","tmt_biascompl","tmt_biashr",usa.ctrls), collapse=" + "), sep=" ~ ")), data = usa.d.comp1and2, Hess = TRUE)

### USA plot for ICC DV, all models (ALL TREATMENTS)
dwplot(list(m1.usa.ols,m1.usa.ols.comp1and2,m1.usa.ctrl.ols, m1.usa.ctrl.ols.comp1and2, m1.usa.log, m1.usa.log.comp1and2, m1.usa.ctrl.log, m1.usa.ctrl.log.comp1and2, m1.usa.olog, m1.usa.olog.comp1and2, m1.usa.ctrl.olog, m1.usa.ctrl.olog.comp1and2),
       vline = geom_vline(xintercept = 0)) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("USA") +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = c(0.007, 0.01),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title.align = .5
    #    ) + ylim(breaks=c("treatment","tmt_bias","tmt_biascompl","tmt_biashr"))  +
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("OLS","OLS Comp1and2","OLS w/ctrl","OLS w/ctrl Comp1and2","Logit","Logit Comp1and2","Logit w/ctrl","Logit w/ctrl Comp1and2","Ord. Logit","Ord. Logit Comp1and2","Ord. Logit w/ctrl","Ord. Logit w/ctrl Comp1and2")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/coeffs_USA_dv1_alltreat_Comp1and2.pdf")


### PHL, ZAF, ISR  (No manip checks in GEO survey)

phl.d$iccsupport_num <- factor(phl.d$iccsupport_num)
zaf.d$iccsupport_num <- factor(zaf.d$iccsupport_num)
isr.d$iccsupport_num <- factor(isr.d$iccsupport_num)

phl.d.comp1and2 <- subset(phl.d, complier1and2 == 1)
zaf.d.comp1and2 <- subset(zaf.d, complier1and2 == 1)
isr.d.comp1and2 <- subset(isr.d, complier1and2 == 1)

### OLS regressions
### No controls, just regression of support on treatment
m1.phl.ols.comp1and2 <-lm(iccsupport_bin ~ treatment, data = phl.d.comp1and2)
m1.zaf.ols.comp1and2 <-lm(iccsupport_bin ~ treatment, data = zaf.d.comp1and2)
m1.isr.ols.comp1and2 <-lm(iccsupport_bin ~ treatment, data = isr.d.comp1and2)

### With "standard" controls, ones that aren't country-specific
m1.phl.ctrl.ols.comp1and2 <-lm(as.formula(paste("iccsupport_bin", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d.comp1and2)
m1.zaf.ctrl.ols.comp1and2 <-lm(as.formula(paste("iccsupport_bin", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d.comp1and2)
m1.isr.ctrl.ols.comp1and2 <-lm(as.formula(paste("iccsupport_bin", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d.comp1and2)

### Logit regressions
### No controls, just regression of support on treatment
m1.phl.log.comp1and2 <-glm(iccsupport_bin ~ treatment, data = phl.d.comp1and2, family = "binomial")
m1.zaf.log.comp1and2 <-glm(iccsupport_bin ~ treatment, data = zaf.d.comp1and2, family = "binomial")
m1.isr.log.comp1and2 <-glm(iccsupport_bin ~ treatment, data = isr.d.comp1and2, family = "binomial")

### With "standard" controls, ones that aren't country-specific
m1.phl.ctrl.log.comp1and2 <-glm(as.formula(paste("iccsupport_bin", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d.comp1and2, family = "binomial")
m1.zaf.ctrl.log.comp1and2 <-glm(as.formula(paste("iccsupport_bin", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d.comp1and2, family = "binomial")
m1.isr.ctrl.log.comp1and2 <-glm(as.formula(paste("iccsupport_bin", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d.comp1and2, family = "binomial")

### Ordered Logit regressions

### No controls, just regression of support on treatment
m1.phl.olog.comp1and2 <-polr(iccsupport_num ~ treatment, data = phl.d.comp1and2, Hess = TRUE)
m1.zaf.olog.comp1and2 <-polr(iccsupport_num ~ treatment, data = zaf.d.comp1and2, Hess = TRUE)
m1.isr.olog.comp1and2 <-polr(iccsupport_num ~ treatment, data = isr.d.comp1and2, Hess = TRUE)

### With "standard" controls, ones that aren't country-specific
m1.phl.ctrl.olog.comp1and2 <-polr(as.formula(paste("iccsupport_num", paste(c("treatment",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d.comp1and2, Hess = TRUE)
m1.zaf.ctrl.olog.comp1and2 <-polr(as.formula(paste("iccsupport_num", paste(c("treatment",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d.comp1and2, Hess = TRUE)
m1.isr.ctrl.olog.comp1and2 <-polr(as.formula(paste("iccsupport_num", paste(c("treatment",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d.comp1and2, Hess = TRUE)

### ISR plot for ICC DV, all models (ALL TREATMENTS)
dwplot(list(m1.isr.ols,m1.isr.ols.comp1and2,m1.isr.ctrl.ols, m1.isr.ctrl.ols.comp1and2, m1.isr.log, m1.isr.log.comp1and2, m1.isr.ctrl.log, m1.isr.ctrl.log.comp1and2, m1.isr.olog, m1.isr.olog.comp1and2, m1.isr.ctrl.olog, m1.isr.ctrl.olog.comp1and2),
       vline = geom_vline(xintercept = 0)) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("ISR") +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = c(0.007, 0.01),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title.align = .5
    #    ) + ylim(breaks=c("treatment","tmt_bias","tmt_biascompl","tmt_biashr"))  +
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("OLS","OLS Comp1and2","OLS w/ctrl","OLS w/ctrl Comp1and2","Logit","Logit Comp1and2","Logit w/ctrl","Logit w/ctrl Comp1and2","Ord. Logit","Ord. Logit Comp1and2","Ord. Logit w/ctrl","Ord. Logit w/ctrl Comp1and2")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/coeffs_ISR_dv1_alltreat_Comp1and2.pdf")

### PHL plot for ICC DV, all models (ALL TREATMENTS)
dwplot(list(m1.phl.ols,m1.phl.ols.comp1and2,m1.phl.ctrl.ols, m1.phl.ctrl.ols.comp1and2, m1.phl.log, m1.phl.log.comp1and2, m1.phl.ctrl.log, m1.phl.ctrl.log.comp1and2, m1.phl.olog, m1.phl.olog.comp1and2, m1.phl.ctrl.olog, m1.phl.ctrl.olog.comp1and2),
       vline = geom_vline(xintercept = 0)) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("PHL") +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = c(0.007, 0.01),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title.align = .5
    #    ) + ylim(breaks=c("treatment","tmt_bias","tmt_biascompl","tmt_biashr"))  +
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("OLS","OLS Comp1and2","OLS w/ctrl","OLS w/ctrl Comp1and2","Logit","Logit Comp1and2","Logit w/ctrl","Logit w/ctrl Comp1and2","Ord. Logit","Ord. Logit Comp1and2","Ord. Logit w/ctrl","Ord. Logit w/ctrl Comp1and2")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/coeffs_PHL_dv1_alltreat_Comp1and2.pdf")

### ZAF plot for ICC DV, all models (ALL TREATMENTS)
dwplot(list(m1.zaf.ols,m1.zaf.ols.comp1and2,m1.zaf.ctrl.ols, m1.zaf.ctrl.ols.comp1and2, m1.zaf.log, m1.zaf.log.comp1and2, m1.zaf.ctrl.log, m1.zaf.ctrl.log.comp1and2, m1.zaf.olog, m1.zaf.olog.comp1and2, m1.zaf.ctrl.olog, m1.zaf.ctrl.olog.comp1and2),
       vline = geom_vline(xintercept = 0)) +
  xlab("Coefficient Estimate") + ylab("") + xlim(-0.7,0.7) +
  ggtitle("ZAF") +
  theme(
    plot.title = element_text(face = "bold"),
    legend.position = c(0.007, 0.01),
    legend.justification = c(0, 0),
    legend.background = element_rect(colour = "grey80"),
    legend.title.align = .5
    #    ) + ylim(breaks=c("treatment","tmt_bias","tmt_biascompl","tmt_biashr"))  +
  ) + ylim(breaks=c("treatment"))  +
  scale_colour_hue(
    name = "Model",
    labels = c("OLS","OLS Comp1and2","OLS w/ctrl","OLS w/ctrl Comp1and2","Logit","Logit Comp1and2","Logit w/ctrl","Logit w/ctrl Comp1and2","Ord. Logit","Ord. Logit Comp1and2","Ord. Logit w/ctrl","Ord. Logit w/ctrl Comp1and2")
  )
ggsave("/Users/stephenchaudoin/Dropbox/Complementarity ZC/coeffs_ZAF_dv1_alltreat_Comp1and2.pdf")










##############################
##################
#### Appendix H ####
##################
##############################

##################
#### Table A9 ####
##################

### Table is just text from the experiment instruments.


##############################
##################
#### Appendix I ####
##################
##############################

##################
#### Table A10 ####
##################

## Regressions 
m1.geo.ols.int <-lm(invsupport_agree ~ treatment + treatment:legalsystem_capable + legalsystem_capable, data = geo.d)
summary(m1.geo.ols.int)

m1.geo.ctrl.ols.int <-lm(as.formula(paste("invsupport_agree", paste(c("treatment + treatment:legalsystem_capable + legalsystem_capable",geo.ctrls), collapse=" + "), sep=" ~ ")), data = geo.d)
summary(m1.geo.ctrl.ols.int)

m1.phl.ols.int <-lm(iccsupport_bin ~ treatment + treatment:legsyscapable_bin + legsyscapable_bin, data = phl.d)
summary(m1.phl.ols.int)

m1.phl.ctrl.ols.int <-lm(as.formula(paste("iccsupport_bin", paste(c("treatment + treatment:legsyscapable_bin + legsyscapable_bin",phl.ctrls), collapse=" + "), sep=" ~ ")), data = phl.d)
summary(m1.phl.ctrl.ols.int)

m1.usa.ols.int <-lm(suppinv_bin ~ treatment + treatment:uslegalsystem_bin + uslegalsystem_bin, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0))
summary(m1.usa.ols.int)

m1.usa.red.ctrl.ols.int <-lm(as.formula(paste("suppinv_bin", paste(c("treatment + treatment:uslegalsystem_bin + uslegalsystem_bin",usa.ctrls), collapse=" + "), sep=" ~ ")), data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0))
summary(m1.usa.red.ctrl.ols.int)

m1.zaf.ols.int <-lm(iccsupport_bin ~ treatment + treatment:otherlegsyscapable_bin + otherlegsyscapable_bin, data = zaf.d)
summary(m1.zaf.ols.int)

m1.zaf.ctrl.ols.int <-lm(as.formula(paste("iccsupport_bin", paste(c("treatment + treatment:otherlegsyscapable_bin + otherlegsyscapable_bin",zaf.ctrls), collapse=" + "), sep=" ~ ")), data = zaf.d)
summary(m1.zaf.ctrl.ols.int)

m1.isr.ols.int <-lm(iccsupport_bin ~ treatment + treatment:isrlegsyscapable_bin + isrlegsyscapable_bin, data = isr.d)
summary(m1.isr.ols.int)

m1.isr.ctrl.ols.int <-lm(as.formula(paste("iccsupport_bin", paste(c("treatment + treatment:isrlegsyscapable_bin + isrlegsyscapable_bin",isr.ctrls), collapse=" + "), sep=" ~ ")), data = isr.d)
summary(m1.isr.ctrl.ols.int)


##################
#### Tables A11 and A12 ####
##################
m1.usa.int <-lm(suppinv_bin ~ treatment + treatment:demrep_ordered + demrep_ordered, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0))
m1.usa.ctrl.int <-lm(suppinv_bin ~ treatment + treatment:demrep_ordered + demrep_ordered + female+age+baorhigher+hhi+region_midwest+region_south+region_west+white_bin++newshours_morethan6, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0))

m1.isr.ols.int <-lm(iccsupport_bin ~ treatment + treatment:leftright + leftright, data = isr.d)
m1.isr.ctrl.ols.int <-lm(iccsupport_bin ~ treatment + treatment:leftright + leftright +female+age+education_postsec+income_26thto75th+income_above75th+jewish+newshours_morethan4, data = isr.d)

m1.phl.ols.int <-lm(iccsupport_bin ~ treatment + treatment:pdpl + pdpl, data = phl.d)
m1.phl.ctrl.ols.int <-lm(iccsupport_bin ~ treatment + treatment:pdpl + pdpl + female + age + education_postsec + income_26thto75th + income_above75th + cebuano + tagalog + catholic + newshours_morethan4 + metromanilancr, data = phl.d)

m1b.phl.ols.int <-lm(iccsupport_bin ~ treatment + treatment:rightparty + rightparty, data = phl.d)
m1b.phl.ctrl.ols.int <-lm(iccsupport_bin ~ treatment + treatment:rightparty + rightparty + female + age + education_postsec + income_26thto75th + income_above75th + cebuano + tagalog + catholic + newshours_morethan4 + metromanilancr, data = phl.d)

m1.zaf.ols.int <-lm(iccsupport_bin ~ treatment + treatment:party_anc + party_anc, data = zaf.d)
m1.zaf.ctrl.ols.int <-lm(iccsupport_bin ~ treatment + treatment:party_anc + party_anc + female + age + education_postsec + income_26thto75th + income_above75th + English + Zulu + Xhosa + christian + party_da + newshours_morethan5 + gauteng, data = zaf.d)

### Tables to match the order in the paper

stargazer(m1.isr.ols.int, m1.isr.ctrl.ols.int, m1.zaf.ols.int, m1.zaf.ctrl.ols.int,		type = "latex",
          title = "Effect of Treatment on Support for ICC Investigation (OLS), w/ Party Interactions, ISR and ZAF",
          label="tab:tmtdomsuppintideolisrzaf",
          keep=c("treatment","leftright","treatment:leftright","treatment:party_anc","party_anc","Constant"),
          dep.var.labels=c("ISR","ISR","ZAF","ZAF"),
          add.lines = c(list(c("Controls?","N","Y","N","Y"))), omit.stat=c("f","ser","adj.rsq","aic","ll"))

stargazer(m1.phl.ols.int, m1.phl.ctrl.ols.int, m1.usa.int, m1.usa.ctrl.int,		type = "latex",
          title = "Effect of Treatment on Support for ICC Investigation (OLS), w/ Party Interactions, PHL and USA",
          label="tab:tmtdomsuppintideolphlusa",
          keep=c("treatment","pdpl","treatment:pdpl","treatment: demrep_ordered","demrep_ordered","Constant"),
          dep.var.labels=c("PHL","PHL","USA","USA"),
          add.lines = c(list(c("Controls?","N","Y","N","Y"))), omit.stat=c("f","ser","adj.rsq","aic","ll"))





##################
#### Table A13
##################
#########################################################
# ICC knowledge summary statistics
#########################################################

# Numerical codings for the icc knowledge variables

isr.d <- isr.d %>% mutate(knowiccbefore.num = case_when(
  knowiccbefore == "A great deal" ~ 5,
  knowiccbefore == "A lot" ~ 4,
  knowiccbefore == "A moderate amount" ~ 3,
  knowiccbefore == "A little" ~ 2,
  knowiccbefore == "None at all" ~ 1
)
)

phl.d <- phl.d %>% mutate(knowiccbefore.num = case_when(
  knowiccbefore == "A great deal" ~ 5,
  knowiccbefore == "A lot" ~ 4,
  knowiccbefore == "A moderate amount" ~ 3,
  knowiccbefore == "A little" ~ 2,
  knowiccbefore == "None at all" ~ 1
)
)

zaf.d <- zaf.d %>% mutate(knowiccbefore.num = case_when(
  knowiccbefore == "A great deal" ~ 5,
  knowiccbefore == "A lot" ~ 4,
  knowiccbefore == "A moderate amount" ~ 3,
  knowiccbefore == "A little" ~ 2,
  knowiccbefore == "None at all" ~ 1
)
)

# USA slightly different question format, and has some NAs "Before today, how much did you know about the ICC?"
usa.d <- usa.d %>% mutate(knowiccbefore.num = case_when(
  q76 == "A lot of information." ~ 4,
  q76 == "Some information." ~ 3,
  q76 == "Very little." ~ 2,
  q76 == "Nothing." ~ 1
)
)


# Binary codings for the icc knowledge variables, 1 if you didn't pick "a little" or "none"

isr.d <- isr.d %>% mutate(knowiccbefore.bin = case_when(
  knowiccbefore.num >= 3 ~ 1,
  knowiccbefore.num < 3 ~ 0
)
)

phl.d <- phl.d %>% mutate(knowiccbefore.bin = case_when(
  knowiccbefore.num >= 3 ~ 1,
  knowiccbefore.num < 3 ~ 0
)
)

zaf.d <- zaf.d %>% mutate(knowiccbefore.bin = case_when(
  knowiccbefore.num >= 3 ~ 1,
  knowiccbefore.num < 3 ~ 0
)
)

usa.d <- usa.d %>% mutate(knowiccbefore.bin = case_when(
  knowiccbefore.num >= 3 ~ 1,
  knowiccbefore.num < 3 ~ 0
)
)


# PHL is the highest self-reported knowledge country, followed by ZAF and ISR.  
library(lattice)
#histogram(phl.d$knowiccbefore.num)
#histogram(zaf.d$knowiccbefore.num)
#histogram(isr.d$knowiccbefore.num)

table(isr.d$knowiccbefore.num)
table(phl.d$knowiccbefore.num)
table(zaf.d$knowiccbefore.num)
table(usa.d$knowiccbefore.num)

summary(isr.d$knowiccbefore.num)
summary(phl.d$knowiccbefore.num)
summary(zaf.d$knowiccbefore.num)
summary(usa.d$knowiccbefore.num)

summary(isr.d$knowiccbefore.bin)
summary(phl.d$knowiccbefore.bin)
summary(zaf.d$knowiccbefore.bin)
summary(usa.d$knowiccbefore.bin)







##################
#### Table A14
##################

# Were uptake rates different in the three countries
manip1.isr <- lm(manip1_correct ~ treatment, data = isr.d)
manip1.phl <- lm(manip1_correct ~ treatment, data = phl.d)
manip1.zaf <- lm(manip1_correct ~ treatment, data = zaf.d)
manip1.usa <- lm(ccheck1_correct ~ treatment, data = usa.d)

summary(manip1.isr)
summary(manip1.phl)
summary(manip1.zaf)
summary(manip1.usa)


manip2.isr <- lm(manip2_correct ~ treatment, data = isr.d)
manip2.phl <- lm(manip2_correct ~ treatment, data = phl.d)
manip2.zaf <- lm(manip2_correct ~ treatment, data = zaf.d)
manip2.usa <- lm(ccheck2_correct ~ treatment, data = usa.d)

summary(manip2.isr)
summary(manip2.phl)
summary(manip2.zaf)
summary(manip2.usa)



##################
#### Table A15
##################

# Interaction terms, does the effect of treatment depend on knowledge and does this differ by country?

usa <-lm(suppinv_bin ~ treatment + treatment:knowiccbefore.bin + knowiccbefore.bin, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0))
usa.ctrl <-lm(suppinv_bin ~ treatment + treatment:knowiccbefore.bin + knowiccbefore.bin + demrep_ordered + female+age+baorhigher+hhi+region_midwest+region_south+region_west+white_bin++newshours_morethan6, data = subset(usa.d, tmt_bias == 0 & tmt_biascompl == 0 & tmt_biashr == 0))
#	USA: knowledge tends to weaken treatment, insignificant

isr <-lm(iccsupport_bin ~ treatment + treatment:knowiccbefore.bin + knowiccbefore.bin, data = isr.d)
isr.ctrl <-lm(iccsupport_bin ~ treatment + treatment:knowiccbefore.bin + knowiccbefore.bin + leftright +female+age+education_postsec+income_26thto75th+income_above75th+jewish+newshours_morethan4, data = isr.d)
#	ISR: treatment _is_ more effective among the more knowledgeable; significant without controls, barely misses 0.1 with

phl <-lm(iccsupport_bin ~ treatment + treatment:knowiccbefore.bin + knowiccbefore.bin, data = phl.d)
phl.ctrl <-lm(iccsupport_bin ~ treatment + treatment:knowiccbefore.bin + knowiccbefore.bin + pdpl + female + age + education_postsec + income_26thto75th + income_above75th + cebuano + tagalog + catholic + newshours_morethan4 + metromanilancr, data = phl.d)
#	PHL: knowledge tends to weaken treatment, insignificant

zaf <-lm(iccsupport_bin ~ treatment + treatment:knowiccbefore.bin + knowiccbefore.bin, data = zaf.d)
zaf.ctrl <-lm(iccsupport_bin ~ treatment + treatment:knowiccbefore.bin + knowiccbefore.bin + party_anc + female + age + education_postsec + income_26thto75th + income_above75th + English + Zulu + Xhosa + christian + party_da + newshours_morethan5 + gauteng, data = zaf.d)
#	ZAF: knowledge tends to make the treatment more negative, interaction term is significant.

### Results table

stargazer(isr, isr.ctrl, zaf, zaf.ctrl, phl, phl.ctrl, usa, usa.ctrl,		type = "latex",
          title = "Effect of Treatment on Support for ICC Investigation (OLS), w/ Prior Knowledge Interactions",
          label="tab:knowinteraction",
          keep=c("treatment","treatment:knowiccbefore.bin","knowiccbefore.bin","Constant"),
          dep.var.labels=c("ISR","ISR","ZAF","ZAF","PHL","PHL","USA","USA"),
          add.lines = c(list(c("Controls?","N","Y","N","Y","N","Y","N","Y"))), omit.stat=c("f","ser","adj.rsq","aic","ll"))



